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Abstract 

Quantum effects due to the spatial delocalization of light atoms are treated in molecular sim¬ 
ulation via the path integral technique. Among several methods, Path Integral (PI) Molecular 
Dynamics (MD) is nowadays a powerful tool to investigate properties induced by spatial delocal¬ 
ization of atoms; however computationally this technique is very demanding. The abovementioned 
limitation implies the restriction of PIMD applications to relatively small systems and short time 
scales. One possible solution to overcome size and time limitation is to introduce PIMD algo¬ 
rithms into the Adaptive Resolution Simulation Scheme (AdResS). AdResS requires a relatively 
small region treated at path integral level and embeds it into a large molecular reservoir consisting 
of generic spherical coarse grained molecules. It was previously shown that the realization of the 
idea above, at a simple level, produced reasonable results for toy systems or simple/test systems 
like liquid parahydrogen. Encouraged by previous results, in this paper we show the simulation of 
liquid water at room conditions where AdResS, in its latest and more accurate Grand-Canonical- 
like version (GC-AdResS), is merged with two of the most relevant PIMD techniques available 
in literature. The comparison of our results with those reported in literature and/or with those 
obtained from full PIMD simulations shows a highly satisfactory agreement. 
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INTRODUCTION 


The structure and dynamics of liquids consisting of molecules that contain light atoms (e.g. 
hydrogen) can be influenced by the quantum effects due to the delocalization of atoms in 
space. In simulation such systems are treated by modeling the atoms of the molecules via 


the path httegta, fototaheat of Feyantaa In patticulat h,utd watet is a typtca, suhiect 

of interest given its role in many helds |^. As explained more in detail in next sections, 
the computational effort is massive because the number of interatomic interactions becomes 
much larger compared to the classical case. As a consequence the size of the system and the 
simulation time affordable with standard computer resources is rather limited. For liquid 
water at room condition a system of 500 molecules for a simulation time of 1-2 ns is usually 
considered already expensive. The limited size and simulation time may imply that particle 
number density fluctuations are arbitrarily suppressed and some systems cannot be treated 
if not at high computational prize (e.g. solvation of a large molecule in water). An opti¬ 
mal complementary technique would consist of a Grand Canonical-like scheme where (local) 
properties can be calculated by employing a computationally affordable path integral simula¬ 
tion of a small open region which, in statistical and thermodynamic equilibrium, exchanges 
particles and energy with a reservoir acting at small computational cost. One possible 
implementation of a Grand Canonical-like Molecular dynamics technique is the Adaptive 


’ ileso 
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ution Simulation scheme (AdResS) 


fl-El 


in its most accurate version of GC-AdResS 


12l |. For the simplest version of AdResS it was shown that for a toy system (liquid of 


tetrahedral molecules) the embedding of a PIMD technique into the scheme produced rather 


encouraging results 


ISj; such results were conhrmed and empowered by the application to 


simple/test systems like liquid parahydrogen at low temperature 15]. In the mean¬ 
while the increased accuracy and more solid conceptual framework of the adaptive scheme 
(GC-AdResS) allows for the study of more complex systems and the calculation of a larger 


number of properties than before 


12l |. In this perspective, this paper reports the technical 


implementation of two different approaches to PIMD, Refs. 


16 


1811 and Refs 


3,Q, 


into 


our GG-AdResS. We show its application to liquid water and report results about static 
and dynamic properties. The comparison with reference data is highly satisfactory and sug¬ 
gest that GC-AdResS, as a complementary method, may play an important role in future 
applications of PIMD (today not feasible with full PIMD simulations). One can think, for 
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example, of solvation of a large molecules (e.g. fullerene in water) and look at possible quan¬ 
tum effects in the structure of the solvation shell. Moreover, GC-AdResS may be employed 
as a tool of analysis and study how the quantum effects change as a function of the size 
of the region treated at PI level. This would represent a novel type of analysis because it 
unequivocally dehnes the essential molecular degrees of freedom required for a given prop¬ 
erty |2l| and thus it allows to quantify how localized are (possible) quantum effects (for the 


properties considered). The paper is organized as follows: next section is dedicated to a 
summary of the relevant technical and conceptual characteristics of GC-AdResS, it follows 
the section dedicated to the description of the basic characteristics of the two PIMD meth¬ 
ods employed in this study. Next the implementation of PIMD in GG-AdResS, for each of 
the two specihc techniques used, is reported. It follows the section of results divided into the 
subsection of (i) static and (ii) dynamic properties. In (i) we report particle number density 
prohles, probability distributions and radial distribution functions of the GG-AdResS simu¬ 
lation compared with results from full PIMD simulations. In (ii) we report the calculation 
of equilibrium time correlation functions compared, also in this case, with data obtained 
from full PIMD simulations. Finally the section of discussion and conclusion is presented. 
The appendix instead reports all technical data of the simulations so that the results can be 
reproduced/checked by other groups. 


GC-ADRESS 

In the original AdResS the coupling idea is rather simple, that is, in a region of interest 
(the atomistic or high resolution region) all the molecular degrees of freedom are treated via 
molecular dynamics while in a (larger) region of minor interest only coarse-grained degrees 
of freedom are treated. The passage of a molecule from one region to another should be 
performed smoothly with a hybrid dynamics in such a way that the atomistic and the coarse¬ 
grained regions are not perturbed in a signihcant way. In order to do so the space is divided 
into three regions, the atomistic (high resolution) region, the coarse-grained region and an 
interfacial region where the atomistic degrees of freedom are transformed in coarse-grained 
and vice versa, we call this region hybrid region or transition region (see FigH]). The coupling 
is made via a space dependent force interpolation: 

= w{X^)w{X^)F^l°^ + [1 - w{X^)w{X^)]F^^ ( 1 ) 
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FIG. 1. Pictorial representation of the GC-AdResS scheme; GG indicates the coarse-grained 
region, HY the hybrid region where path-integral and coarse-grained forces are interpolated via a 
space-dependent, slowly varying, function w{x) and EX (or PI) is the path-integral region (that 
is the region of interest). Top, the standard set up with the thermostat that acts globally on 
the whole system, used in the calculation of static properties. Bottom, the “local” thermostat 
technique employed in this work in the calculation of dynamical properties. 


where a and jS indicate two molecules, and w{Xa) and w{Xis) indicate the interpolating 
(weighting) function depending on the coordinate of the center of mass of the molecules X^ 
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and Xp\ 


X < (Iat 


w(x) = < 


cos 




d-AT ^ X <i cIat d/\ 

0 d-AT + dA < X 

where, dAx and dA are size of atomistic and hybrid regions respectively . is the force 

in the atomistic region, which is derived from atomistic interactions, is the force in the 
coarse-grained region, which is derived from a coarse-grained potential. A thermostat takes 
care of thermally eqnilibrating the atomistic degrees of freedom reintroduced in the transition 
region. This simple set up turned out to be computationally robust; the calculation of 
structural and thermodynamics property in AdResS compared with the calculations done 


in a subregion of equivalent size in_a 
agreement for several test systems 


21 


nil atomistic simulation shows a highly satisfactory 


26| . The computational robustness encouraged the 
the method on the basis of hrst principles of 
This analysis hrst led to the introduction 


investigation of the conceptual justihcation o 
thermodynamics and statistical mechanics 
of a thermodynamic force acting on the center of mass of the molecules in the hybrid region. 
The thermodynamic force is based on the principle of uniformizing, to the atomistic value, 
the chemical potential of each (space dependent) resolution and then to the derivation 
of such a thermodynamic force from a more general thermodynamic principle, that is from 
the balance of grand potential for two interfaced open systems Q : 


Patom T Po I d^th{,X^dT 


V = PcgV 


( 2 ) 


where Patom and Pcg are the pressure of the atomistic and coarse-grained region, po is the 
target density of the reference full atomistic simulation, V the volume of the simulation box. 
The explicit calculation of Fth{x) is reported in the next section. Based on such derivation a 
step forward was done and AdResS was reformulated in terms of Grand Canonical formalism 
(GC-AdResS) where mathematical rigorous conditions were derived in order to assure that 
in the atomistic region the system samples a Grand Canonical distribution. Such condition, 
at the hrst order, have been shown to be equivalent to the use of the thermodynamic force 
[lO, 11|. Moreover the coarse-grained model can be arbitrarily chosen, without any reference 


to the atomistic model. Recent results 


11 


Ensemble model of Bergmann and Lebowitz 


rave embedded the scheme into the Grand 


29 


30[ | and introduced a local thermostat acting 


only in the coarse-grained and hybrid region. Such a formalization allows one to dehne, with 
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well founded physical arguments, the Hamiltonian of the atomistic (high resolution) region 
as the kinetic energy plus the interaction energy of the molecules in the atomistic region 
only; this implies that the interaction with molecules outside can be formally neglected. The 
dehnition of the Hamiltonian allows then to properly dehne the procedure for the calculation 
of equilibrium time correlation functions; moreover, for the case of PI approach, this set up 
will provide a rigorous dehnition of the Hamiltonian of quantization. As it will be specihed 
later on, there exists also a clear numerical argument that supports the dehnition of an 
accurate Hamiltonian in the high resolution region. 


PIMD TECHNIQUES 


The path integral formalism of Feynman applied to molecular simulation/dynamics of molec¬ 
ular systems is a well established approach and thus here we will not report its formal deriva¬ 
tion but only those aspects which are technically relevant for this specihc study. A formal 


derivation and discussion of basic aspect of this approach can be found, in Refs.[16|,lMi for ex¬ 


ample. The essential point of interest (in this paper) is the transformation, via path integral 
formalism, of a classical Hamiltonian of N distinguishable particles with phase space coor¬ 
dinate (p, r), mass rrij (for the j-th particle) and interaction potential in space U(ri, ....r^v): 

N 2 

= E (3) 

j=l J 

into a quantized Hamiltonian which is formally equivalent to a Hamiltonian of classical 
polymer rings (atoms). The interatomic potential is distributed over the beads in such a 
way that each bead of a polymer ring interacts with the corresponding bead of another 
polymer ring. The intra-atomic interactions consists of harmonic potentials which couple 
each bead to the first neighbors in the chain. The fictitious dynamics of this polymeric 
liquid, with the spatial fluctuations/oscillations of the rings describing the quantum spatial 
delocalization of the atoms, allows for the calculation of quantum statistical properties of 
the atomic/molecular system. The quantized Hamiltonian takes the form: 


i=l 


N 


ip'*')? 


N 




• 1 

J = 1 J 


2 / (i) 


i=i 


‘"F + pV(A, -■i-t) 


(4) 


where P is the number of beads of the polymer, and p* are a fictitious mass 

VP 

/3n 


and momentum respectively, up = ^ (/? = l/kpT) and V{vp ....r^y) is the potential that 
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acts between same bead index i of two different particles. This set up allows to use molec¬ 
ular dynamics for the calculation of statistical properties. However the direct use of the 
Hamiltonian above has shown to lead to a highly non-ergodic dynamics and suffers from 
poor sampling problems in the extended phase space of polymer ring [^, since there are a 
wide range of frequencies present. The highest frequency limits the time step to be used in 
the simulation which causes the low frequency modes to be poorly sampled. Thus either a 
very small time step or very long runs should be performed, starting from different initial 
conditions in order to overcome this problem. In order to circumvent the ergodicity prob¬ 
lem, normal modes transformation is preferred 3, The basic idea is to decouple the 
harmonic spring term, so that only a single harmonic frequency remains in the dynamics, 
and the time step for the simulation can be adjusted accordingly. The whole procedure is 
based on a transformation of coordinates to normal mode coordinates and thus to the use 
of an effective Hamiltonian: 


Hp = 


i=l 


■ N {if 


N 




where is the potential that acts between same bead index i of 


j=i 2m 
dh/ 


h) 


Si) I 


(5) 


two different particles in terms of the normal mode coordinates .Xjy. 


A. Choice of masses 


In the standard PIMD 


33 


34l |. the masses are chosen such that all the internal modes 


have the same frequency and the sampling is efficient. Thus the choice of mass is: 

.P 


=m,A!-,i = 2, 


= mj 


where mj is the physical mass and A® are the eigenvalues obtained by the normal mode 
transformation. This approach was used to calculate static properties and here we will use 
it, within GC-AdResS, for the same purpose. We will refer to this approach as HI. Craig 


and Manolopolous 


35| have developed ring polymer molecular dynamics (RPMD), which 


has been successfully shown to calculate time correlation functions; the choices of the masses 
in RPMD is as follows: 

( 6 ) 


(>) 

mj = mj 
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In this work, we will employ this approach within GC-AdResS to calculate, in addition to 


static properties, time correlation functions; we wil 
there exists an alternative formulation for RPMD 


refer to it as H2 approach. However, 


19| . The classical Hamiltonian for RPMD 


is: 


2 = 1 


N 


y lP_ii + y 

^ 2m, ^ 




N 


rrii 


,(*) 


Ai - r, 


P+1)n2 


r + nrl,....r]v) 


(7) 


where (3p = fi/P., which effectively means that the simulation is performed at P times 
the original temperature. Moreover the harmonic bead-bead interaction and the potential 
energy are scaled by P relative to EqO In Ref. 18|, equivalence between different RPMD 
formalisms was shown. Due to the calculation of the thermodynamic force, for GC-AdResS 
simulations this becomes an interesting technical aspect to investigate (see next sections). 
We will refer to this approach as H3 and verify its numerical robustness in GG-AdResS by 
comparison with the results obtained from H1,H2. 


B. PIMD in GC-AdResS 


The original idea of merging PIMD and AdResS was based on a simple extension of the 
AdResS principle. The dynamics of polymer rings, from a technical point of view, is nothing 
else than the dynamics of classical degrees of freedom, thus the standard AdResS could be 
applied (technically) in the same way, with only the modification 13l-ll5|: 


= w{X^)w{X^)F^^ + [1 - w{X^)w{X^)]F^^ 


( 8 ) 


where is the force between beads of the rings representing the atoms of molecule a 
and molecule 13. However one of the authors has shown before that in any adaptive scheme, 
based on a spatial interpolation of atomistic and coarse grained interactions, it cannot exist a 


valid rigorous global Hamiltonian 


36| . Thus from the conceptual point of view the coupling 


between the polymer rings and the coarse-grained molecules cannot be rigorously expressed 
in a Hamiltonian form. However, calculations have shown that PIMD-AdResS was able to 
reproduce very well results obtained with full PIMD simulations. Since the Hamiltonian 


formalism is at the basis of the PIMD approach the procedure of Refs. 13l-ll5| was empirical 


and could be verified only a posteriori. The reason why the procedure was successful is 
that the coupling between the polymer rings and the coarse-grained molecules is negligible, 
in terms of energetic contribution, under the hypothesis that the path integral region and 






















the coarse-grained region were large enough compared to the hybrid region. However we 
have also numerically verified that even when all the three regions are relatively small and 
comparable in size, results are still satisfactory. The latest formalization of AdResS in GC- 
AdResS, reported in the previous section, justifies why from a conceptual point of view the 
setup of PI-AdREsS is robust. In fact according to the model of Bergmann and Lebowitz 




for a simulation in a Grand Ensemble one does not need to have an explicit 
coupling between the path integral region and the reservoir. The necessary and sufficient 
condition is the knowledge of the molecules’ distribution in the reservoir. It follows that 
the interaction of the molecules of the path integral region with the rest of the system, 
while technically convenient and numerically efficient, from the conceptual (formal) point of 
view instead does not play a crucial role. Such interaction plays only the technical role of 
a sort of “capping potential” which avoids that molecules entering the path integral region 
overlap in space. Moreover the action of the thermodynamic force and of the thermostat in 
the hybrid region makes the stochastic coupling dominant (compared to the explicit hybrid 
interactions), which is the essence of any Grand Ensemble scheme. It follows that in GG- 
AdResS the Hamiltonian to be considered for the path integral formalism is the Hamiltonian 
of the path integral region only, without any external additional term; i.e. the path integral 
region with its quantized Hamiltonian is embedded in a large reservoir with the proper 
Grand Ganonical behaviour. It must be clarified that while the Bergmann-Lebowitz model 
provides an elegant and solid formal structure to the PI-AdResS, however it is not strictly 
required to justify the existence of an accurate Hamiltionian in the PI region and thus the 
implementation of PIMD in AdREsS. In fact in the appendix we provide a numerical proof 
that, for the systems treated in this paper, the interaction energy between the PI region 
and the rest of the system is at least one order of magnitude smaller than the interaction 
energy of the molecules in the PI region. The accuracy and robustness of PI-AdResS (or 
PI-GG-AdREsS) will be shown with the simulation of liquid water in the next section. 
Finally it must be clarified that for the current implementation of PIMD in GG-AdResS (in 
the GROMAGS package), it is difficult to estimate the computational gain since the code 
architecture is not yet optimized. At this stage we want only to show that the approach is 
satisfactory from a conceptual point of view. However, for very large systems with P = 32, 
the computational gain is around 1.7-2.0 compared to the full PIMD simulations. With 
further code modifications (e.g. removal of explicit degrees of freedom in the coarse-grained 
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region, using multiple time steps) or with the implementation of PI-AdResS in a platform 
explicitly designed for PIMD simulation we estimated, for systems of the order of thousand 
molecules, a gain of at least a factor 4.0-5.0 compared to the full PIMD simulations. 


1. Calculation of the Thermodynamic Force in PIMD 


For an atomistic system, the thermodynamic force, Fth{x), can be expressed as: 

Fth{x) = —VP{x) (9) 

Po 

where M is the mass of the molecule and P{x) is the pressure which characterizes each 
different resolutions (for the initial conhguration). P{x) is approximated in terms of linear 
interpolation of molecular number density: 

M 

P{x) = Patom H- [Po - p(a:)] (10) 

poKi 

where p{x) is the density generated if the simulation runs without any thermodynamic force. 
The thermodynamic force is then obtained by an iterative procedure: 

- 7^vp,(i) (11) 

[poV^^ 

After each iteration, a density prohle p{x) is obtained due to the application of the ther¬ 
modynamic force. The process converges when the density prohle obtained is equal to the 
target density. At this point the system is in thermodynamic equilibrium and the production 
run can start. The calculation of thermodynamic force in PIMD-GC-AdResS is essentially 
based on the same principle of balancing grand potential for interfaced open systems: 


Pquantum T Po / Fth{r)dr 


V = PcgV 


( 12 ) 


L JA J 

where po is the target density of the reference full path-integral system. As for the classical 
case, P{x) can be written as: 


R(x) Pquantum T [Po P(^)] 


(13) 


While the above approach is highly efficient for classical simulations, for path integral sim¬ 
ulations, it is cumbersome to run an PIMD-GC-AdResS simulations to calculate the ther¬ 
modynamic force, before doing an actual production run, as the path-integral simulation is 
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inherently very expensive. In order to make the scheme efficient we have devised a strategy 
to calculate the thermodynamic force which requires least computation. As discussed in the 
previous section, we will show how the thermodynamic force is calculated for the different 
Hamiltonian approaches. In case of HI, and H2 where the temperature of the system is just 
the normal temperature, we calculated the thermodynamic force for path-integral systems 
with varying Trotter number P = 1,4, 6,8 and 10 (P = 1 represents the classical limit). 
Since the thermodynamic force takes care of a thermodynamic equilibration and since the 
thermodynamic conditions (thermodynamic state point) of a classical and a quantum sys¬ 
tem are the same, we expect that the thermodynamic force calculated in the classical case 
(P = 1) is sufficient to provide thermodynamic equilibrium in simulations where P = 32 
is used. In fact we found that the thermodynamic force was same in all the cases. Fig [2] 
shows the thermodynamic force calculated for a water system, with different number of ring 
polymer beads in each case. Using this argument, we used this thermodynamic force in the 
actual production run with P = 32. We found that the density of water molecules in the 
full quantum subregion and the transition region is equal to the reference density of the 
water system at the same thermodynamic conditions. Thus, in the HI and H2 approach, 
if the quantum effects on the pressure of the system are not large, we can directly use the 
thermodynamic force calculated from the classical simulation. 


In H3 approach, the situation is more complex, as the effective temperature of simulation 
changes if the number of beads is changed, thus the (numerical) thermodynamic state point 
changes. In this case, there would be no other choice but to run a full PIMD-GC-AdResS 
simulation with P = 32 and calculate the thermodynamic force. However, we avoided 
such an expensive calculation and instead calculated the thermodynamic force for system 
with different number of beads P = 1,4, 6, 8 and 10 at temperatures T = 298 x P and 
extrapolated thermodynamic force for P = 32, using space dependent factors calculated 
from thermodynamic force for smaller values of P. Next we used this thermodynamic force 
for production run with P = 32, and found that the density of water molecules in the full-PI 
subregion is same as the target density while the density in transition region differs at worst 
by 3%. 
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FIG. 2. Thermodynamic force calculated in AdResS simulation using HI (H2) approach. The 
force is calculated for different number of polymer ring beads. It does not change as the number 
of beads is varied. 

2. Equilibrium Time Correlation Functions: Theoretical and computational aspects 


The technique of Ring Polymer molecular dynamics (RPMD) (H2) focuses on the Kubo- 


transformed correlation functions 




the operators A and B is dehned by 


1 


35j: 


The Kubo-transformed correlation function of 




where Z is the canonical partition function: 


Z = tr 




(14) 


(15) 


The RPMD method approximates the Kubo-transformed correlation functions by using the 
classical ring-polymer trajectories generated by the dynamics produced by the Hamiltonian 
in Eq. [71 The RPMD approximation is given by 39| : 


CAsit) 


{2'KhfP^Zt 
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d’’Pad’’Y, A‘Ara)B'Ar,) 


(16) 


2=1 
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FIG. 3. Thermodynamic force calculated in AdResS simulation using H3 approach. The force is 
calculated for different number of polymer ring beads. The thermodynamic force for P = 32 is 
then extrapolated by using space-dependent scaling factors calculated using thermodynamic force 
for P = 1,4,6, 8 and 10. 

where Zp is the canonical partition function, and r* indicates the time evolution at time t 
of the positions. The functions Ap{ro) and Bp{rt) are calculated by taking the average over 
the beads of the ring polymer: 

P P 

Mr) = pY. M),Bp(r) = i S(r) (17) 

i=i i=i 

For the calculations in GC-AdResS the above equation needs to be written in the formalism 
of the Grand Ganonical ensemble: 

, " ( 18 ) 
^ J^J2Mr«mBUr,{r»{N))) 

i=l 

where p is the chemical potential and N' is the number of molecules at time ‘O’, that 
remain correlated at time ‘t’ (that is molecules which remain in the path integral region 
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for the whole time within the time window considered); Zp^ = is the grand- 

canonical partition function and Hp{N) is the Hamiltonian of the (open) path integral 
region with N, instantaneous number of molecules. It must be noticed that the a priori 
knowledge of pi is not required, actually in GC-AdResS p, is automatically calculated by the 
equilibration procedure of the thermodynamic force (see also [h]). From the technical point 
of view we have used the same calculation procedure as that of Ref. , where equilibrium 
time correlation functions were calculated in the open subsystems using classical molecular 
dynamics. Such a principle is based on the dehnition of reservoir in the Bergmann-Lebowitz 
model, which implies that when a molecule leaves the system and enters the reservoir, it 
looses its microscopic identity and thus the corresponding correlations; thus, if a molecule 
which is present at time to, disappears from the system at time t (i.e. moves into the 
reservoir), then the contribution of this molecule, outside the time window [0,f], to the 
correlation function shall not be considered. In our previous work we have shown that 
such a principle is physically consistent on the basis of results of molecular simulations. 
Since all the beads in a ring-polymer are treated as dynamical variables [^, there are 
no thermostats used in RPMD simulations. Thus, the simulations are performed under 
NVE conditions, with either starting conhgurations generated from massively-thermostated 
PIMD simulations [^, or re-sampling of momenta from Maxwell-Boltzmann distribution 


after every few picoseconds 4l|]. In order to keep the dynamics of the beads Newtonian m 
the path-integral subregion of GC-AdResS, we use a “local-thermostat” procedure HQ, 
where the thermostat is applied only in the coarse-grained and hybrid region, while the 
explicit path-integral region is thermostat-free. This ensures, that the molecules which are 
present in the path-integral subregion are not subject to any perturbation due to the action 
of the thermostat. 


RESULTS 


In this section we report results about the simulation of liquid water at room conditions. 

n 

The quantum model for liquid water used in this work is q-SPC/FW [^. It was shown 
that the thermodynamic and dynamical properties calculated using this water model agree 
quite well with the experiment data. The section is divided in two parts, the hrst where 
static results (molecular number density across the system, radial distribution functions. 
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probability distribution of the molecules) are reported, and the second where several equi¬ 
librium time correlation functions are calculated. Few further points must be mentioned as 
clarification to this study. The total volume of the PIMD-GC-AdResS box is the same in 
all simulations, while three different sizes of the region at PI resolution are used and the 
dimension of the transition region is kept always the same. The smallest size of the PI region 
represents the limiting case of a statistically relevant number of molecules treated with PI 
resolution. The largest size instead represents the limiting case of a reservoir (hybrid plus 
coarse-grained region) which is relatively small and thus it may be expected to not fulhll 
the conceptual requirement of being statistically large enough. We will show that even in 
these two limiting cases the method is computationally and conceptually robust. A second 
point to take into account is that we compare the results of GC-AdResS for the PI region 
with the results obtained in a subsystem of a full PI simulations, such a subsystem is of the 
same size of the GG-AdResS simulation. The subsystem of a large full PI simulation box 
is a natural Grand Ganonical ensemble, thus if our subsystem of AdResS reproduces the 
results of a full PI subsystem then we can be rather confident that the PI region in AdResS 
sample the Grand Ganonical distribution sufficiently well. From the physical point of view, 
it should be clarified that the functions calculated in a subsystem must be considered local 
in space and time if compared to calculation done over the whole simulation box of the full 
PI simulation. Once again, as the subsystem size increases the functions go to the value 


obtained in a fu’ 
checks in Ref. 


1 PI simulation when the full box is considered (physical consistency, see 


12|). Technical details of the simulation are reported in the Appendix. 


C. Static Properties 


We use the HI and H2 PIMD approaches (Hl-GG-AdResS and H2-GG-AdResS respec¬ 
tively for the GG-AdResS simulation), FigJU shows molecular number density. In all three 
cases the agreement is highly satisfactory, the largest deviation is found for the case with PI 
region of 0.5nm and is below 5%. This is the basic test to show equilibration and thermo¬ 
dynamic consistency, moreover, following the mathematical formalization of Ref. lOj] is the 
hrst order necessary condition in order to have the correct Grand Ganonical distribution in 
the PI region. A further conhrmation of the fact that the method samples the phase space 
of a subsystem in a sufficiently correct way is represented by Fig 13 The figure shows the 
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FIG. 4. Molecular number density calculated with GC-AdResS for different size of quantum sub- 
region. Results are compared with the density obtained in a full path integral simulation. 


particle number probability distribution in quantum subregion of AdResS and an equivalent 
subregion in full path integral simulation. It can be seen that also in this case the results 
are highly satisfactory and the shape of two curves is a Gaussian, as one should expect. 
The g(r) is an important structural quantity that represents a two-body correlation function 
and thus a higher order than the molecular density of the ensemble many-body distribution; 
moreover it differs considerably when quantum models of water are used, in particular corre¬ 


lation functions involving hydrogen atoms 


43|. We calculated the local bead-bead g(r)’s 
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FIG. 5. Particle number probability distribution of GC-AdResS compared with the equivalent full 
path integral subsystem, for different size of quantum subregion. The shape of both curves is a 
Gaussian (reference black continuous curve) in all the three different simulations. The top part of 
the figure indicates the extension of the PI region (compared to the rest of the system) where the 
function is calculated; this representation is equivalent in all subsequent figures. 

the quantum subregion in GC-AdResS and compared them with the bead-bead g(r)’s in an 
equivalent subregion in the full path-integral simulation. Fig [B] to Fig [S] show that the results 
from GC-AdResS agree with the results from full PI simulation in a highly satisfactory way. 
We have also verihed, for the most relevant case [EX = 1.2), that also the the H3 approach 
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gives satisfactory results for the static properties when employed in GC-AdResS; results are 
reported in Fig O Due to the more empirical calculation of the thermodynamic force in 
H3 results are not as accurate as those of HI and H2; the density in the hybrid region 
differs by around 3%, which is anyway numerically negligible (however the difference must 
be reported). However, the number probability distribution and bead-bead g(r)’s agree quite 
well in AdResS and full path-integral simulations. This leads to the conclusion that also 
results obtained with H3 are highly satisfactory. 

This section essentially show the ability of PI-GC-AdResS with all three PIMD techniques 
to sample basic (but highly relevant) static properties of a Grand Canonical ensemble. In 
order to prove that a more elaborated sampling is also satisfactorily made by the method 
we report in the next section the calculation of equilibrium time correlation functions. 


D. Dynamic properties 


We report results for the velocity-velocity autocorrelation function, for the hrst and second 


order orientational (molecular dipole) correlation 
correlation function for hydrogen bond dynamics 


unction 
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45| and for the reactive flux 


4y, l47|. This latter, in specihc situations 


may strongly diverge from the classical case, and thus it may be a quantity of relevance; 
moreover the fact that PI-GC-AdResS reproduce the behaviour of a full PI simulation is of 
high technical relevance in perspective (e.g. study of solvation of molecules). The explicit 


formulas used for the functions calculated here are given in Ref. 


5l| . All results shown 


m 


this section are highly satisfactory, either when H2 is used or H3 is used. Thus the PI- 
GC-AdResS can be certainly considered a robust computational method for the calculation 
of quantum-based static and dynamic properties of liquid water and as a consequence for 
simpler systems and for system where water play a major role (at least). 


1. Equilibrium Time Correlation Functions 

Figure [IHl E] and m show the three correlation functions calculated in the quantum sub- 
region in GC-AdResS and an equivalent subregion in RPMD simulation, where the explicit 
region is 1.2 nm. All the correlation functions are calculated using H2 approach and H3, 
which conhrm the consistency of the two methods in GC-AdResS. As stated before, these 
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FIG. 6. From left to right: (bead-bead) oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen 
partial radial distribution functions calculated with path integral AdResS. Such functions are 
compared with the results obtained for an equivalent subsystem (EX = 0.5nm) in a full path 
integral simulation. 


are the local time correlation functions, calculated in the specihc region of interest, and 


could differ from the global time correlation functions, calculated over the whole system. 
However, it was shown in Ref 12| , that as the size of the explicit region increases, the local 


correlation functions converge to the global correlation functions. 
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FIG. 7. From left to right: (bead-bead) oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen 
partial radial distribution functions calculated with path integral AdResS. Such functions are 
compared with the results obtained for an equivalent subsystem (EX = 1.2nm) in a full path 
integral simulation. 


2. Dynamics of hydrogen bonding 


In order to investigate the dynamics of hydrogen bond formation and breaking using RPMD 


simulations we calculate the hydrogen bond population fluctuations in time, which are char- 
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FIG. 8. From left to right: (bead-bead) oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen 
partial radial distribution functions calculated with path integral AdResS. Such functions are 
compared with the results obtained for an equivalent subsystem (EX = 2Anm) in a full path 
integral simulation. 


acterized by the correlation function: 


c(t) = {h(Q)h(t)) / (h) 


(19) 
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FIG. 9. From left to right (top): Particle number probability distribution of GC-AdResS obtained 
using the H3 approach. 

From left to right (bottom): (bead-bead) oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen 
partial radial distribution functions calculated with path integral AdResS using the H3 approach. 
Such functions are compared with the results obtained for an equivalent subsystem (EX = 1.2nm) 
in a full path integral simulation. 
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where h{t) is the hydrogen bond population operator, which has a value 1, when a particular 
pair are bonded, and zero otherwise. One can then calculate the rate of relaxation as: 

k{t) = —dc/dt (20) 

k[t) is the average rate of change of hydrogen-bond population for those trajectories where 
the bond is broken at a time t later. The two water molecules are treated as hydrogen bonded, 
if the distance between the center of two oxygen rings is less than 0.35 nm and simultaneously 
the angle between the axis dehned by the center of two oxygen ring polymers, and center of 
one of the oxygen-hydrogen ring is less than SOdegrees. Fig [13] shows k{t) calculated in the 
quantum subregion of AdResS and an equivalent subregion in RPMD simulation. 



FIG. 10. Kubo-transformed velocity auto correlation function for q-SPC/FW water model calcu¬ 
lated in quantum subregion of GC-AdResS and an equivalent subregion in RPMD simulation. 


CONCLUSION 

We have performed simulations of liquid water at room conditions using PIMD in three 
different technical approaches. Each of these approaches was embedded in GC-AdResS so 
that a PIMD for open systems in contact with a generic reservoir is realized. The results 
regarding static and dynamic quantities is highly satisfactory and qualihed PI-GC-AdResS 
as a robust method for simulations of systems which currently are prohibitive for full PIMD 
simulations. For example the already mentioned solvation problem. One can dehne a high 
resolution region at PI resolution around the solute and surround the solvation region with 
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FIG. 11. Kubo-transformed first order orientational correlation function for q-SPC/FW water 
model calculated in quantum subregion of GC-AdResS and an equivalent subregion in RPMD 
simulation. Dipole moment axis is chosen as the inertial axis of molecule. 



FIG. 12. Kubo-transformed second order orientational correlation function or q-SPC/FW water 
model calculated in quantum subregion of GC-AdResS and an equivalent subregion in RPMD 
simulation. Dipole moment axis is chosen as the inertial axis of molecule. 

a reservoir as that constructed in GC-AdResS. The static and dynamic properties of the 
hydrogen bonding network can be analyzed and, by comparing results with those of classical 
systems, one may conclude about the importance of quantum effects due to hydrogen spatial 
delocalization. This approach can introduce not only a technical innovation regarding the 
computational efficiency, but, by varying the size of the high resolution region, could be 
used as a tool of analysis to identify the essential degrees of freedom required by a certain 
physical process. In this perspective, here we have shown that PI-GC-AdResS is a robust 
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FIG. 13. The rate function k(t) for q-SPC/FW water model calculated in quantum subregion of 
GC-AdResS and an equivalent subregion in RPMD simulation. 


method for linking the microscopic to macroscopic scale in a truly multiscale fashion. 
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I. APPENDIX: TECHNICAL DETAILS 


A. Energetic contribution of the coupling term 


The i-th molecule (at position, Uj) in the EX (PI) region is characterized by w{ri) = 1. It 
follows that the force acting on the Tth molecule can be separated in two parts; (i) the force 
generated by the interaction of molecule i with molecules of the EX region: 

F,j = F^J,Wj e EX (21) 

and (ii) the force generated by the interaction with molecules in the rest of the system: 

F,,, = f/ + |1 - w(ri)]Fff,Vj € HY + CG. (22) 

From EqJ2I]it follows: 

(23) 

where V, is the gradient w.r.t. molecule i and Upj is a compact form to indicate the proper 
bead-bead interaction of atoms of molecule i with those of molecule j. Eql22] represents 
instead the coupling force between molecules of HY + CG region and molecule i, that is an 
external force. At this point we argue that the non-integrable part of the dynamics in the 
HY region is a numerically negligible boundary effect. In fact Eql22]can be rewritten as: 

E + E + [1 - <i’(>-,)lv,c/cGl. 

jeHY+CG j&HY+CG 

(24) 

It follows that the energy of the i-th molecule at a certain time time t associated with the 
force of EqlH] is given by: 


j&HY+GG 

where the Res = HY -|- CG. The total energy of coupling at time t is then defined as: 

wp,.a,.m = E ■ (26) 

i&pi 

In order to understand whether or not the quantity of Eqj26] is numerically negligible, one 
should compare it to the amount of energy, Wpi^pj, corresponding to the interaction between 
molecules of the PI region only: Wpi=pj{t) = Ylii<j Upi'CC ^ PI- If 

\Wpj-Pl{t) \ — \ Wpj-Res{t)\ 


|hPp7-P/(t)| 


l;Vt 


(27) 
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FIG. 14. Main figure: Wpi-pi(t) compared to Wpi-Res(t). Inset: The relative amount 
of the interaction between the PI region and the rest of the system along the trajectory: 

contribution is, at most, of 10%. Calculations are done within the 
HI and H2 approach (top) and H3 (bottom). 

then it seems reasonable to approximate the total energy of the PI region by the Hamiltonian 
of the PI region, thus the Hamiltonian formalism is numerically justihed in PI-AdResS. 
FigHH shows that the difference in energy is at least of one order of magnitude and that 
condition [27] holds in all simulations we have presented in this work. Moreover it should be 
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noticed that on purpose we have performed simulations where the technical conditions are 
not optimal (the size of each region of the system is much smaller than the size prescribed 
by the theory), thus Eql27] would certainly hold in simulations with standard technical 
conditions. 


B. simulation set up 


Static Properties: All path integral simulations are performed by home-modified GRO- 
MACS 1491, and the thermodynamic force in GC-AdResS simulations is calculated using 
VOTGA j^. The number of water molecules in system are 1320, and the box dimensions 
are 5.8 x 2.6 x 2.6 corresponding to a density 990 kgm~^. In AdResS simulations, 
the resolution of the molecules changes along x-axis, as depicted in Figure [H Three dif¬ 
ferent AdResS simulations are performed, each with a different size of quantum subregion. 
The different sizes of the quantum subregion treated in this work are 0.5 x 2.6 x 2.6 nrrfi, 
1.2 X 2.6 X 2.6 nm? and 1.2 x 2.6 x 2.6 nm?. The transition region, which has dimen¬ 
sions 2.6 X 2.6 X 2.6 nrrfi is fixed in all the three cases. The remaining system contains 
coarse-grained particles, which interact via generic WGA potential of the form: 


U{r) = 4e 


12 


6n 


+ e,r < 2^/® 


cr 


( 28 ) 


The parameters a and e in the current simulations are 0.30 nm and 0.65 kJ/mol respectively. 
Thirty two ring polymer beads are used in all the simulations, which is sufficient to obtain 
converged results for both static and dynamical properties. Reaction field method is used 
to compute the electrostatic properties with dielectric constant for water = 80. The cut-off 
for both van der Waals and electrostatic interactions is 1.2 nm. All the static properties 
are computed from 250 ps long trajectories. The simulations using HI and H2 formalisms 
are performed at 298 iP, while the simulations using iJ3 formalism are performed at 9536 
K. The time step used in all the simulations is 0.1 fs. In the calculation of thermodynamic 
force, a single iteration consists of a 200 ps long trajectory which is used to compute the 
density profile. A total of 20 such iterations is sufficient to obtain a fiat density profile, and 
a converged thermodynamic force. 

Dynamic Properties: The system details are kept same as in the previous section. A 200 
ps long PIMD simulation is performed and along the trajectory, configurations are taken 






after every 8 ps to perform RPMD simulations. Thus a total of 25 trajectories each of 
length 25 ps are generated. For the hrst 5 ps, we keep the thermostat switched on, in order 
to adjust the velocities as masses are different in PIMD and RPMD methods. After this 
initial equilibration run, the thermostat is switched off, and the NVE trajectories generated 
are used to compute various time correlation functions. We use the same strategy for 
AdResS simulations, where a 200 ps long fully thermostated GC-AdResS PIMD simulation 
is performed, and 25 initial conhgurations are taken along this trajectory to perform GC- 
AdResS RPMD simulations. For the hrst 5 ps, the thermostat acts in the explicit as well as 
the hybrid and coarse-grained regions. After the short equilibration run, the thermostat is 
switched off in the explicit region, while the hybrid and coarse-grained region are kept under 
the action of the thermostat. The dynamic properties are calculated in the explicit region in 
the last 20 ps, i.e. excluding the equilibration run. The velocity auto-correlation function is 
calculated for 1 ps, while the orientational correlation functions and reactive hux correlation 
functions for hydrogen bond dynamics are calculated for 10 ps in one single trajectory, and 
then averaged over all the trajectories. 


C. thermostat issue 


It is well known that massive thermosetting is needed in the path integral simulations, as 
the forces arising due to the high frequencies in the polymer ring and the forces due to 
the potential U{x) are weakly coupled. Tuckerman et ah coupled each normal mode 
variable to separate Nose-Hoover chains, thereby ensuring proper ergodic sampling of the 
phase space. Manolopolous et. al [19| developed specihc Langevin equations for thermostat, 
that are tuned to sample all the internal modes of the ring polymer quite efficiently. However, 
in this work we chose the standard Langevin equations of thermostat with time scale 0.1 
ps, which is strong enough for sampling the phase space effectively, though it may not be 
be the most efficient choice. The reason is that in the initial stage of validating GC AdResS 
for path integral simulations, we need to show that the properties obtained in the full PI 
simulations are reproduced exactly in AdResS. Since we use the same thermostat in both the 
simulations, there should not be any discrepancy arising due to the thermostat. However, 
the comparison of static properties calculated in our reference PIMD simulation with those 
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available in literature (referring to the approaches above) is highly satisfactory. 
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